R Notebook

Il primo step sarà quello di installare tutti i pacchetti necessari:

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("curatedTCGAData")
BiocManager::install("TCGAutils")
BiocManager::install("TCGAbiolinks")
install.packages("SNFtool")
install.packages("NetPreProc")
install.packages("caret");
install.packages("cluster");
install.packages("mclustcomp");

quindi carico le librerie:

library("curatedTCGAData");
library("TCGAbiolinks");
library("TCGAutils");
library("SNFtool");
library("NetPreProc");
library("caret");
library("cluster"); #pam
library("mclustcomp");

A questo punto per ottenere il dataset del carcinoma prostatico (PRAD) contenente informazioni su espressione genica a livello di mRNA, espressione di microRNA e espressione proteica, ho seguito il seguente processo:

# ho specificato le omiche di interesse richieste dal progetto
assays <- c("miRNASeqGene", "RNASeq2Gene", "RPPAArray");

# ottengo i dati per il dataset del carcinoma prostatico
mo <- curatedTCGAData(diseaseCode = "PRAD", 
                        assays = assays, 
                        version = "2.0.1", dry.run = FALSE);

snapshotDate(): 2023-10-24
Working on: PRAD_miRNASeqGene-20160128
see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation
loading from cache
Working on: PRAD_RNASeq2Gene-20160128
see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation
loading from cache
Working on: PRAD_RPPAArray-20160128
see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation
loading from cache
Working on: PRAD_colData-20160128
see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation
loading from cache
Working on: PRAD_metadata-20160128
see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation
loading from cache
Working on: PRAD_sampleMap-20160128
see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation
loading from cache
harmonizing input:
  removing 5189 sampleMap rows not in names(experiments)

mo

A MultiAssayExperiment object of 3 listed
 experiments with user-defined names and respective classes.
 Containing an ExperimentList class object of length 3:
 [1] PRAD_miRNASeqGene-20160128: SummarizedExperiment with 1046 rows and 547 columns
 [2] PRAD_RNASeq2Gene-20160128: SummarizedExperiment with 20501 rows and 550 columns
 [3] PRAD_RPPAArray-20160128: SummarizedExperiment with 195 rows and 352 columns
Functionality:
 experiments() - obtain the ExperimentList instance
 colData() - the primary/phenotype DataFrame
 sampleMap() - the sample coordination DataFrame
 `$`, `[`, `[[` - extract colData columns, subset, or experiment
 *Format() - convert into a long or wide DataFrame
 assays() - convert ExperimentList to a SimpleList of matrices
 exportClass() - save data to flat files

Fase 2: Data pre-processing Ogni campione è identificato da un barcode con una specifica struttura: - i primi 12 caratteri identificano uno specifico individuo - le altre parti ci danno indicazioni sul tipo di campione (primario, metastatico, solido, derivato dal sangue, ecc.), sul tipo di materiale genomico estratto (DNA, RNA) e altre informazioni relative alle repliche tecniche (cioè misurazioni ripetute dallo stesso campione).

Usiamo il barcode per: - conservare solo i tumori solidi primari per avere un gruppo di campioni più omogeno. Questo è identificato dal codice “01” nella parte “sample” del barcode - Rimuovere eventuali duplicati rappresentati da campioni che hanno gli stessi 12 caratteri iniziali

# Considerare solo i tumori solidi primari:
primary <- TCGAutils::TCGAsampleSelect(colnames(mo), c("01"))

Avvertimento: Inconsistent barcode lengths: 28, 27Avvertimento: Inconsistent barcode lengths: 28, 27

mo <- mo[, primary, ]

harmonizing input:
  removing 106 sampleMap rows with 'colname' not in colnames of experiments

# Controllo delle repliche 
check_rep <- anyReplicated(mo)
print(check_rep)

PRAD_miRNASeqGene-20160128  PRAD_RNASeq2Gene-20160128    PRAD_RPPAArray-20160128 
                     FALSE                      FALSE                      FALSE 

# Verifica se ci sono valori nulli in ciascuna colonna di un dataframe
for (i in 1:length(complete)){
  print(paste(names(complete)[i], ":", any(colSums(is.na(complete[[i]])) > 0)))
}

[1] "PRAD_miRNASeqGene-20160128 : FALSE"
[1] "PRAD_RNASeq2Gene-20160128 : FALSE"
[1] "PRAD_RPPAArray-20160128 : FALSE"

# Remove features having NAs (present only in proteomics data):
complete[[3]] <- complete[[3]][, colSums(is.na(complete[[3]])) == 0];

# Remove features with near zero variance and retain top 500 features 
# having higher variance:
nf <- 100;
for(i in 1:length(complete)){
    
    idx <- caret::nearZeroVar(complete[[i]])
    message(paste("Removed ", length(idx), "features from", names(complete)[i]));
    if(length(idx) != 0){
        complete[[i]] <- complete[[i]][, -idx];
    }

    if(ncol(complete[[i]]) <= nf) next
    
    vars <- apply(complete[[i]], 2, var);
    idx <- sort(vars, index.return=TRUE, decreasing = TRUE)$ix;
    
    complete[[i]] <- complete[[i]][, idx[1:nf]];
    as.data.frame(complete[[i]])
    
}

Removed  0 features from PRAD_miRNASeqGene-20160128
Removed  1334 features from PRAD_RNASeq2Gene-20160128
Removed  0 features from PRAD_RPPAArray-20160128

# Perform features standardization using z-score:
zscore <- function(data){
    
    zscore_vec <- function(x) { return ((x - mean(x)) / sd(x))}
    data <- apply(data, 2, zscore_vec)
    
    
    return(data)
}

complete <- lapply(complete, zscore);

# Clean barcodes retaining only "Project-TSS-Participant":
for(v in 1:length(complete)){
    rownames(complete[[v]]) <- substr(rownames(complete[[v]]), 1, 12);
}

PUNTO 3:

# Print number of samples for each subtype:
table(subtypes$Subtype_Integrative);

  1   2   3 
 60  83 105 

# Verifica che i pazienti nei dataset multi-omics e nelle sottocategorie siano nello stesso ordine
matching_order <- all(!is.na(match(rownames(complete[[1]]), substr(subtypes$pan.samplesID, 1, 12))))

if (matching_order) {
  cat("I pazienti nei dataset multi-omics e nelle sottocategorie sono nello stesso ordine.\n")
} else {
  cat("ATTENZIONE: I pazienti nei dataset multi-omics e nelle sottocategorie NON sono nello stesso ordine.\n")
}

I pazienti nei dataset multi-omics e nelle sottocategorie sono nello stesso ordine.

PUNTO 4:

matching_order <- TRUE

for (i in 1:length(complete)) {
  matching_order <- matching_order & !is.na(match(rownames(complete[[i]]), substr(subtypes$pan.samplesID, 1, 12)))
}

if (all(matching_order)) {
  cat("I pazienti nei dataset multi-omics e nelle sottocategorie sono nello stesso ordine dopo il riordinamento.\n")
} else {
  cat("ATTENZIONE: I pazienti nei dataset multi-omics e nelle sottocategorie NON sono nello stesso ordine.\n")
}

I pazienti nei dataset multi-omics e nelle sottocategorie sono nello stesso ordine dopo il riordinamento.

K = 20; # Number of neighbors
T = 20; # Number of iterations

# Compute similarity matrix for each data source using the scaled
# exponential euclidean distance:
W_list <- list();
for(i in 1:length(complete)){
    Dist <- (dist2(as.matrix(complete[[i]]), as.matrix(complete[[i]])))^(1/2);
    W_list[[i]] <- affinityMatrix(Dist);
}
    
# Integration of multi-omics data using Similarity Network Fusion:
# t is the number of iterations and K is the number of neighbors to 
# consider to compute the local similarity matrix:
W_int <- SNF(W_list, K, T)

PUNTO 6:

# Calcolo della media delle matrici di similarità
W_average <- Reduce(`+`, W_list) / length(W_list)

PUNTO 8a:

# b. Normalizza le matrici di similarità
for (i in 1:length(W_list)) {
    W_list[[i]] <- NetPreProc::Max.Min.norm(W_list[[i]])
}

# c. Converti le matrici di similarità normalizzate in matrici di distanza
D_list_single <- lapply(W_list, function(W) 1 - NetPreProc::Max.Min.norm(W))

# Esegui PAM sugli insiemi di dati integrati
k <- length(unique(subtypes$Subtype_Selected))

# PAM per miRNA
pam_res_miRNA <- pam(as.dist(D_list_single[[1]]), k = k)

# PAM per mRNA
pam_res_mRNA <- pam(as.dist(D_list_single[[2]]), k = k)

# PAM per proteine
pam_res_proteins <- pam(as.dist(D_list_single[[3]]), k = k)

PUNTO 8b:

# Normalizza la matrice di similarità
W_avg_norm <- NetPreProc::Max.Min.norm(W_average)

# Calcola la matrice di distanza
D_avg <- as.dist(1 - W_avg_norm)

# Esegui PAM sulla matrice integrata ottenuta dalla media
pam_res_avg <- pam(D_avg, k = k)

PUNTO 8C:

# Normalizza la matrice di similarità fusionata
W_int_norm <- NetPreProc::Max.Min.norm(W_int)

# Calcola la matrice di distanza
D_int <- as.dist(1 - W_int_norm)

# Applica PAM sull'insieme di dati integrato tramite Similarity Network Fusion
pam_res_snf <- pam(D_int, k = k)

PUNTO 10 (Opzionale):

spectral_clustering <- spectralClustering(W_int, k)
str(spectral_clustering)

 int [1:248] 1 1 3 1 3 3 2 2 2 2 ...

# Visualizza i risultati
table(spectral_clustering)

spectral_clustering
 1  2  3 
83 84 81 

clustering_results <- list(pam_res_miRNA = pam_res_miRNA$clustering,
                           pam_res_mRNA = pam_res_mRNA$clustering,
                           pam_res_proteins = pam_res_proteins$clustering,
                           pam_res_avg = pam_res_avg$clustering,
                           pam_res_snf = pam_res_snf$clustering,
                           spectral_clustering = spectral_clustering)

# Calcolare la tabella di contingenza per ciascun risultato di clustering rispetto ai sottotipi di malattie iCluster
for (name in names(clustering_results)) {
  print(paste("Contingency table for", name))
  print(table(subtypes$Subtype_Integrative, clustering_results[[name]]))
}

[1] "Contingency table for pam_res_miRNA"
   
     1  2  3
  1 24 15 21
  2 51 23  9
  3 59 30 16
[1] "Contingency table for pam_res_mRNA"
   
     1  2  3
  1 15 18 27
  2  7 26 50
  3 42 18 45
[1] "Contingency table for pam_res_proteins"
   
     1  2  3
  1 24 23 13
  2 38 32 13
  3 33 35 37
[1] "Contingency table for pam_res_avg"
   
     1  2  3
  1 29 25  6
  2 37 24 22
  3 25 44 36
[1] "Contingency table for pam_res_snf"
   
     1  2  3
  1 14  9 37
  2 50 18 15
  3 14 65 26
[1] "Contingency table for spectral_clustering"
   
     1  2  3
  1 15 16 29
  2 12 28 43
  3 56 40  9

# Calcolare l'indice di Rand Adjusted per ciascun risultato di clustering rispetto ai sottotipi di malattie iCluster
for (name in names(clustering_results)) {
  print(paste("Adjusted Rand index for", name))
  print(cluster::adjustedRandIndex(iCluster_subtypes, clustering_results[[name]]))
}

[1] "Adjusted Rand index for pam_res_miRNA"

Errore: 'adjustedRandIndex' non è un oggetto esportato dal 'namespace:cluster'

PUNTO 11:

# Creare una tabella di confronto con nomi descrittivi per gli approcci
comparison_table <- data.frame(
  Method = c("PAM_miRNA", "PAM_mRNA", "PAM_proteins", "PAM_Avg", "PAM_SNF", "Spectral_Clustering"),
  Cluster_Assignment = c(
    pam_res_miRNA$clustering, 
    pam_res_mRNA$clustering, 
    pam_res_proteins$clustering, 
    pam_res_avg$clustering, 
    pam_res_snf$clustering, 
    spectral_clustering
  )
)

# Visualizzare la tabella
print(comparison_table)